clear all;
% clc;

% ------------------------- Geometry Parameters ---------------------------

% import data
patientgeom = load('femspine_patientgeom.txt');
cordgeom = load('patient_cordgeom.txt');

patnum = 5;
ptag=['p',num2str(patnum),'_'];
noaxons = 200;

epidur  = 0; % 0 = intradural, 1 = epidural
subloc  = 1; % 1 = 1 mm above cord; 2 = 1 mm below dura
fval = [cordgeom(patnum,8),cordgeom(patnum,9),1-cordgeom(patnum,9)];
tval = [0,-10,-20];
nofv  = length(fval);
notv  = length(tval);
% polarity of stimulation
pol = -1;
if(pol>=0)
    poltag = 'pos_';
else
    poltag = 'neg_';
end

epi_tag = [ptag,poltag,'epi_bp_dc_act'];
sub_tag = [ptag,poltag,'sub_bp_dc_act'];

epi_drtag = [ptag,poltag,'epi_bp_dr_act'];
sub_drtag = [ptag,poltag,'sub_bp_dr_act']; 

% ------------------------ Import Epidural Data ---------------------------

epibp_data  = zeros(noaxons,3);  
epidr_data  = zeros(noaxons,3);
s1 = strrep(num2str( fval(1) ),'.','p');
cd(['patient',num2str(patnum),'_thresholds']);

for jj = 1:notv
    
    s2 = strrep(num2str(tval(jj)),'-','n');
    loctag = ['_d',s1,'_t',s2];
      
    file_jj   = [epi_tag,loctag,'.txt'];
    drfile_jj = [epi_drtag,loctag,'.txt'];
    epibp_temp = load(file_jj);
    epidr_temp = load(drfile_jj);
    
    epibp_data(:,jj) = sort(epibp_temp(:,1)*pol);
    epidr_data(:,jj) = sort(epidr_temp(:,1)*pol);
    
end

% ------------------------ Import Subdural Data ---------------------------

subbp_data = zeros(noaxons,6);
subdr_data = zeros(noaxons,6);

for ii = 1:2
    
    s1 = strrep( num2str( fval(ii+1) ),'.','p' );
 
    for jj = 1:notv
        
        kk = notv*(ii-1)+jj;
        s2 = strrep( num2str( tval(jj) ),'-','n' );     
        loctag = ['_d',s1,'_t',s2];
        subbp_temp = load([sub_tag,loctag,'.txt']);
        subdr_temp = load([sub_drtag,loctag,'.txt']);
        
        subbp_data(:,kk) = sort(subbp_temp(:,1)*pol);
        subdr_data(:,kk) = sort(subdr_temp(:,1)*pol);
        
    end
    
end
cd ..;

% ------------------------ Importing Currents -----------------------------

epi_iunit = zeros(1,3);
sub_iunit = zeros(1,6);

cd(['patient',num2str(patnum),'_potentials']);
s1 = strrep(num2str( fval(1) ),'.','p');
for jj = 1:notv
    
    s2 = strrep(num2str(tval(jj)),'-','n');
    loctag = ['_d',s1,'_t',s2];
      
    file_jj = [ptag,'epi_bp_current',loctag,'.txt'];
    epi_iunit(jj) = load(file_jj);

end

for ii = 1:2
    
    s1 = strrep( num2str( fval(ii+1) ),'.','p' );
 
    for jj = 1:notv
        
        kk = notv*(ii-1)+jj;
        s2 = strrep( num2str( tval(jj) ),'-','n' );     
        loctag = ['_d',s1,'_t',s2];
        
        file_kk = [ptag,'sub_bp_current',loctag,'.txt'];
        sub_iunit(kk) = load(file_kk);
        
    end
    
end
cd ..;

% ------------------------ Processing Data --------------------------------

xp     = 1;
indx_p = round(noaxons/xp);

epi_I = (epi_iunit'*ones(1,noaxons))';
sub_I = (sub_iunit'*ones(1,noaxons))';

epibp_min = min(epibp_data.^2.*(1e3*epi_I),[],2);
epibp_max = max(epibp_data.^2.*(1e3*epi_I),[],2);
tmp1 = epidr_data.^2.*(1e3*epi_I);
epidr_mu10 = mean( tmp1(indx_p,:) );

subbp_min = min(subbp_data.^2.*(1e3*sub_I),[],2);
subbp_max = max(subbp_data.^2.*(1e3*sub_I),[],2);
tmp2 = subdr_data.^2.*(1e3*sub_I);
subdr_mu10 = mean( tmp2(indx_p,:) );

epibp_pow=epibp_data.^2.*(1e3*epi_I);
subbp_pow=subbp_data.^2.*(1e3*sub_I);
relp50=median(subbp_max)/median(epibp_min)-1;
relp50

% ------------------------ Plotting Data ----------------------------------

ylv = 0:1:100;

fact = 100*(1/noaxons:1/noaxons:1);
fact = fact';

figure;
hold on;
c1 = 0.1;
c2 = 0.55;
fill(cat(1,epibp_min,flipud(epibp_max)),cat(1,fact,flipud(fact)),[c1,c1,c1],'LineWidth',1);
fill(cat(1,subbp_min,flipud(subbp_max)),cat(1,fact,flipud(fact)),[c2,c2,c2],'LineWidth',1);
% plot(epidr_mu10*ones(size(ylv)),ylv,'k-','LineWidth',4);
% plot(subdr_mu10*ones(size(ylv)),ylv,'k--','LineWidth',4);
hold off;
% ttl = ['Bipolar Stimulation with AD-TECH Array (Patient ',num2str(patnum),')'];
% title(ttl,'FontSize',30,'FontWeight','b');
xlabel('Stimulation Power (mW)','FontSize',36);
ylabel('DC Fibers Activated (%)','FontSize',36);
L1 = 'Epidural';
L2 = 'Intradural';
L3 = '10 % DR (Epidural)';
L4 = '10 % DR (Intradural)';
%legend(L1,L2,L3,L4,'Location','NW');
legend(L1,L2,'Location','NW');
set(gca,'FontSize',36,'XScale','log','XTick',[0.01,0.1,1,10,100],'YTick',25:25:100);
xlim([0.01,300]);
ylim([0,100]);
axis square;